summary(d)
d$ArrivalTime.mean     <- ave(d$ArrivalTime, d$Nest)
d$ArrivalTime.centered <- d$ArrivalTime - d$ArrivalTime.mean
d$NegPerChick.log10    <- log10(d$NegPerChick + 1)
ggpairs(d[,c(2, 3, 9, 10)]) +
theme_few()
ggplot(data = d, aes(x = ArrivalTime.centered, y = NegPerChick.log10)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~Nest) +
theme_few()
ggpairs(d[,c(2, 3, 8, 9, 10)]) +
theme_few()
beyondoptimal.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = Owls)
library(lme4)
beyondoptimal.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = Owls)
beyondoptimal.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
beyondoptimal.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
plot(beyondoptimal.lmer)
library(predictmeans)
beyondoptimal.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
residplot(beyondoptimal.lmer)
beyondoptimal.glmer <- glmer(NegPerChick ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d, family = poisson)
plot(beyondoptimal.glmer)
residplot(beyondoptimal.glmer)
residplot(beyondoptimal.lmer)
randomintercept.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(1 | Nest), data = d)
fixedeffects.lm <- lm(LogNeg ~ (SexParent + FoodTreatment + mean.arrival) * arrival.centered, data = Owls)
randomintercept.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(1 | Nest), data = d)
fixedeffects.lm <- lm(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3, data = d)
REMLRatioTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
return((pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2)
}
REMLRatioTest(randomslope.lmer, randomintercept.lmer, df = 2) # df = 2 for variance component of slope and its covariance with intercept
randomslope.lmer <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
REMLRatioTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
return((pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2)
}
REMLRatioTest(randomslope.lmer, randomintercept.lmer, df = 2) # df = 2 for variance component of slope and its covariance with intercept
REMLRatioTest(randomintercept.lmer, fixedeffects.lm, df = 1)  # df = 1 for variance component of intercept
randomslope.lmer.ml <- lmer(NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d, REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
library(MuMIn)
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
plot(dd)
bestmodel.lmer <- lmer(NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered +
(ArrivalTime.centered | Nest), data = d)
residplot(bestmodel.lmer)
summary(bestmodel.lmer)
confint(bestmodel.lmer)
?confint
?Confint
library(car)
?Confint
confint(bestmodel.lmer, oldNames = FALSE)
bestmodel.lmer <- lmer(NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered +
(ArrivalTime.centered || Nest), data = d)
residplot(bestmodel.lmer)
bestmodel.lmer <- lmer(NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered +
(ArrivalTime.centered | Nest), data = d)
residplot(bestmodel.lmer)
r.squaredGLMM(best.lmer)
r.squaredGLMM(bestmodel.lmer)
library(effects)
plot(effect("FoodTreatment:ArrivalTime", bestmodel.lmer), multiline = TRUE)
Anova(bestmodel.lmer)
plot(effect("FoodTreatment:ArrivalTime.centered", bestmodel.lmer), multiline = TRUE)
plot(allEffects(bestmodel.lmer))
ranef(best.lmer)
ranef(bestmodel.lmer)
plot(allEffects(bestmodel.lmer))
plot(effect("FoodTreatment:ArrivalTime.centered", bestmodel.lmer), multiline = TRUE)
plot(allEffects(bestmodel.lmer))
plot(effect("FoodTreatment:ArrivalTime.centered", bestmodel.lmer), multiline = TRUE, confint = TRUE)
plot(allEffects(bestmodel.lmer))
plot(allEffects(abund.bestmodel),
main = NULL,
xlab = "Arrival time (centered by nest)",
ylab = expression(log[10](Negotiations per chick)))
plot(allEffects(abund.bestmodel),
main = NULL,
xlab = "Arrival time (centered by nest)",
ylab = expression(log[10](Negotiationsperchick)))
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Arrival time (centered by nest)",
ylab = expression(log[10](Negotiationsperchick)))
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Arrival time (centered by nest)",
ylab = expression(log[10]("Negotiationsperchick")))
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Arrival time (centered by nest)",
ylab = expression(log[10]("Negotiations per chick")))
?S
10^-1.684e-01
S(bestmodel.lmer, brief = TRUE)
?S
0.412 + residplot(bestmodel.lmer, group = "FoodTreatment")
plot(stdres(bestmodel.lmer) ~ d$FoodTreatment)
stdres(bestmodel.lmer)
resid(bestmodel.lmer)
?residuals.merMod
plot(resid(bestmodel.lmer, scaled=TRUE) ~ d$FoodTreatment)
10^0.168
10^-0.168
confint(bestmodel.lmer, oldNames = FALSE)
10^-0.207196073
10^-0.12942820
effects(bestmodel.lmer)
allEffects(bestmodel.lmer)
?allEffects
effect("FoodTreatment", bestmodel.lmer)
10^effect("FoodTreatment", bestmodel.lmer)
effect("FoodTreatment", bestmodel.lmer)
10^0.4115296
10^0.2431531
1.750464/2.579465
effect("ArrivalTime.centered", bestmodel.lmer)
?effect
allEffects(bestmodel.lmer)
??ggeffects
d <- read.delim(file = "Data/owls.txt")
d$Nest <- factor(d$Nest)
d$FoodTreatment <- factor(d$FoodTreatment)
d$SexParent <- factor(d$SexParent)
summary(d)
d <- read.delim(file = "Data/owls.txt")
d$Nest <- factor(d$Nest)
d$FoodTreatment <- factor(d$FoodTreatment)
d$SexParent <- factor(d$SexParent)
d$ArrivalTime.mean     <- ave(d$ArrivalTime, d$Nest)
d$ArrivalTime.centered <- d$ArrivalTime - d$ArrivalTime.mean
d$NegPerChick.log10    <- log10(d$NegPerChick + 1)
ggplot(data = d, aes(x = ArrivalTime.centered, y = NegPerChick.log10)) +
facet_wrap(~Nest) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
library(ggplot2)
library(ggthemes)
ggplot(data = d, aes(x = ArrivalTime.centered, y = NegPerChick.log10)) +
facet_wrap(~Nest) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
library(lme4)
randomslope.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
residplot(randomslope.lmer)
library(predictmeans)
randomslope.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(1 | Nest), data = d)
fixedeffects.lm <- lm(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3,
data = d)
randomslope.lmer.ml <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered || Nest), data = d, REML = FALSE, na.action = "na.fail")
library(MuMIn)
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
plot(dd)
bestmodel.lmer <- lmer(
NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered + SexParent + (ArrivalTime.centered | Nest),
data = d)
S(bestmodel.lmer)
library(car)
S(bestmodel.lmer)
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^2 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(1 | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
S(g1)
S(randomslope.glmer )
randomslope.glmer2 <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(1 | Nest), data = d, family = "poisson")
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
REMLBoundaryTest(randomslope.glmer1, randomslope.glmer2, 2)
library(ggeffects)
?ggpredict
# Inspect the data
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$Sex <- factor(d$Sex)
summary(d)
ggpairs(d) +
theme_few()
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
ggpairs(d) +
theme_few()
ggpairs(d[,-3]) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex + (age | Subject), data = d)
residplot(randomslope.lmer)
ggpairs(d[,-3]) +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex + age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age, data = d)
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, na.action = "na.omit", REML = FALSE)
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.omit")
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.omit")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
plot(dd)
residplot(randomintercept.lmer.ml)
# Re-fit the reduced model (using REML)
bestmodel.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
# Check assumptions of the reduced model
residplot(bestmodel.lmer)
# Interpret the reduced model
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
ggpairs(d[,-3]) +
theme_few()
ggpairs(d[,-3]) +
theme_few()
# Plot distance as a function of age with a panel for each subject, coloured by sex
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
residplot(randomslope.lmer)
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
ggpairs(d[,-3]) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age, data = d)REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
# Test the significance of the random intercept
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
# Model comparison for fixed effects (using ML)
# We will now re-fit the random-intercept model using ML and dredge the fixed effects
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d,
REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
plot(dd)
residplot(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
library(emmeans)
allEffects(bestmodel.lmer)
effect("age", bestmodel.lmer)
effect(Sex, bestmodel.lmer)
effect("Sex", bestmodel.lmer)
Anova(bestmodel.lmer, test = "F", type =3)
1.816^2
1.386^2
3.297856 + 1.920996
1.816^2 / (1.816^2 + 1.386^2)
S(bestmodel.lmer)
1.816^2 / (1.816^2 + 1.386^2)
1.816^2 / (1.816^2 + 1.386^2)
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$FoodTreatment <- factor(d$FoodTreatment)
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$FoodTreatment <- factor(d$FoodTreatment)
d$Sex <- factor(d$Sex)
summary(d)
# Center age
d$age.mean     <- ave(d$age, d$Subject)
d$age.centered <- d$age - d$age.mean
ggpairs(d[,-(2:3)]) +
theme_few()
# Tutorial 10 Answers
# Load packages
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
ggpairs(d[,-(2:3)]) +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
S(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
S(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
randomintercept.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age.centered, data = d)
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
randomintercept.lmer.ml <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d,
REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
bestmodel.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
residplot(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
Anova(bestmodel.lmer, test = "F", type = 3)
bestmodel.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
Anova(bestmodel.lmer, test = "F", type = 3)
S(bestmodel.lmer)
bestmodel.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
S(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
4.795e-01+3.048e-01
0.2964574+0.0669912
library(emmeans)
emtrends(model, ~ Sex, var = "age.centered")
library(emmeans)
emtrends(bestmodel.lmer, ~ Sex, var = "age.centered")
emmeans(bestmodel.lmer)
emmeans(bestmodel.lmer, ~Sex)
residuals.mermode
library(lme4)
?residuals
??resid
library(predictmeans)
?residplot
log10(1)
?residplot
# Tutorial 10 Answers
# Load packages
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
# Load the data, convert to factors, and inspect the data
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$Sex <- factor(d$Sex)
summary(d)
# Center age
d$age.mean     <- ave(d$age, d$Subject)
d$age.centered <- d$age - d$age.mean
# Plot the data as a scatterplot matrix (excluding the Subject column)
ggpairs(d[,-(2:3)]) +
theme_few()
# Plot distance as a function of age with a panel for each subject, coloured by sex
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
# Build the full model
# We will build a random-slope model including fixed effects of sex, age, and their interaction
# and random effects comprising a random slope and intercept for each subject
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
# Check assumptions of the full model
residplot(randomslope.lmer, group = "Sex")
residplot(randomslope.lmer, group = "Sex", slope = TRUE)
randomslope.resid <- residuals(randomslope.lmer, scaled = TRUE)
plot(randomslope.resid ~ Sex, data = d)
plot(randomslope.resid ~ age.centered, data = d)
plot(randomslope.resid ~ Subject, data = d)
